These analyses investigate connectivity dynamics of the posterior medial network (PMN) during a movie-watching task. The data includes 136 younger adults (age 18-40) (68 for the initial ‘discovery’ analyses, 68 for replication) from the CamCan dataset.

Each PMN region is a cluster of 100 contiguous voxels (2x2x2mm) within the Default A and C subnetworks (Schaefer et al., 2018) that activate during ‘episodic’ tasks (according to the NeuroSynth meta analysis) - based on the definition in Ritchey & Cooper (2020, TiCS).

The analyses first test time-averaged functional connectivity, so the pearson correlation (and partial correlation) between ROI time series during the full movie watching task. I use this to investigate the presence of subsystems (from seed-to-voxel analyses), the most dominant connections within the PMN, and the most influential regions.

Next, the analyses test time-varying connectivity of the PMN subsystems in relation to the movie content (similarity across subjects), including change in connectivity at event transitions.

Setup

library(tidyverse)
library(plotly)
library(ggnetwork)
library(cowplot)
library(ggpubr)
library(knitr)
library(psych)
library(rstatix)
library(reshape2)
library(pracma)
library(mosaic)
library(NetworkToolbox)
library(pander)
library(lessR)
library(fmri)
library(ppcor)
library(corpcor)
library(network)
library(GGally)
library(igraph)
library(tidygraph)
library(lsa)
library(prodlim)
library(stringr)
library(caret)


## define sample to analyze (discovery or replication)
this.group <- "Discovery"
cat(this.group,'sample')
## Discovery sample

R package information

sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] caret_6.0-86         prodlim_2019.11.13   lsa_0.73.1          
##  [4] SnowballC_0.6.0      tidygraph_1.1.2      igraph_1.2.2        
##  [7] GGally_1.4.0         network_1.15         corpcor_1.6.9       
## [10] ppcor_1.1            MASS_7.3-51          fmri_1.9.3          
## [13] lessR_3.8.1          pander_0.6.3         NetworkToolbox_1.2.1
## [16] mosaic_1.4.0         Matrix_1.2-14        mosaicData_0.17.0   
## [19] ggformula_0.9.0      ggstance_0.3.1       lattice_0.20-35     
## [22] pracma_2.1.8         reshape2_1.4.3       rstatix_0.6.0       
## [25] psych_1.8.10         knitr_1.20           ggpubr_0.1.8        
## [28] magrittr_1.5         cowplot_0.9.3        ggnetwork_0.5.8     
## [31] plotly_4.8.0         forcats_0.4.0        stringr_1.4.0       
## [34] dplyr_0.8.5          purrr_0.3.2          readr_1.1.1         
## [37] tidyr_1.0.2          tibble_3.0.2         ggplot2_3.1.0       
## [40] tidyverse_1.2.1     
## 
## loaded via a namespace (and not attached):
##  [1] readxl_1.1.0         backports_1.1.2      plyr_1.8.4          
##  [4] lazyeval_0.2.2       splines_3.5.1        digest_0.6.18       
##  [7] foreach_1.4.4        htmltools_0.3.6      wesanderson_0.3.6   
## [10] cluster_2.0.7-1      mosaicCore_0.6.0     openxlsx_4.1.0      
## [13] recipes_0.1.13       modelr_0.1.2         gower_0.2.1         
## [16] aws_2.4-1            colorspace_1.4-1     rvest_0.3.2         
## [19] ggrepel_0.8.0        haven_1.1.2          crayon_1.3.4        
## [22] jsonlite_1.5         survival_3.1-11      iterators_1.0.10    
## [25] glue_1.3.2           gtable_0.3.0         ipred_0.9-9         
## [28] V8_1.5               car_3.0-2            DEoptimR_1.0-8      
## [31] abind_1.4-5          scales_1.0.0         randomcoloR_1.1.0   
## [34] Rcpp_1.0.1           viridisLite_0.3.0    foreign_0.8-71      
## [37] awsMethods_1.1-1     stats4_3.5.1         lava_1.6.7          
## [40] htmlwidgets_1.3      httr_1.3.1           RColorBrewer_1.1-2  
## [43] ellipsis_0.3.0       pkgconfig_2.0.2      reshape_0.8.8       
## [46] nnet_7.3-12          tidyselect_1.0.0     rlang_0.4.5         
## [49] munsell_0.5.0        cellranger_1.1.0     tools_3.5.1         
## [52] cli_1.1.0            generics_0.0.2       sas7bdat_0.5        
## [55] broom_0.7.0          evaluate_0.12        ggdendro_0.1-20     
## [58] yaml_2.2.0           ModelMetrics_1.2.2.2 zip_1.0.0           
## [61] robustbase_0.93-3    nlme_3.1-137         leaps_3.0           
## [64] xml2_1.2.0           compiler_3.5.1       rstudioapi_0.8      
## [67] curl_4.3             png_0.1-7            stringi_1.4.3       
## [70] gsl_1.9-10.3         vctrs_0.2.4          pillar_1.4.4        
## [73] lifecycle_0.2.0      data.table_1.11.8    R6_2.4.0            
## [76] latticeExtra_0.6-28  gridExtra_2.3        rio_0.5.10          
## [79] codetools_0.2-15     assertthat_0.2.1     rprojroot_1.3-2     
## [82] withr_2.1.2          metafor_2.4-0        mnormt_1.5-5        
## [85] triangle_0.11        parallel_3.5.1       hms_0.4.2           
## [88] grid_3.5.1           rpart_4.1-13         timeDate_3043.102   
## [91] class_7.3-14         rmarkdown_1.10       carData_3.0-2       
## [94] Rtsne_0.15           pROC_1.16.2          lubridate_1.7.4     
## [97] ellipse_0.4.1

Define custom functions & color palettes

### define functions:
se <- function(x) sqrt(var(x)/length(x))  #function to calculate SE
ci <- function(x) (sqrt(var(x)/length(x))) * 1.96  #function to calculate 95% CI

### color palettes:
# color my 8 nodes
mycolors = c('#f2a68c','#934be0','#929292','#e6664d','#8ce9f7','#3c7ff5','#fcd841','#95df3a')
module.colors <- c("#fcd841","#3c7ff5","#8ce9f7")  #for "Core", "Ventral", "Ventral-to-Core"

# color scales
gray.colors  <- colorRampPalette(c("gray90","gray10"))
heat.colors  <- colorRampPalette(c("#009FFF","white","#ff5858"))
blue.colors  <- colorRampPalette(c("gray70","gray90","#56CCF2","#2F80ED","#0575E6"))
phase.colors <- c(heat.colors(13)[13],heat.colors(13)[8],heat.colors(13)[6],heat.colors(13)[1])

Load in event boundaries

# create event boundary windows
my.TR = 2.47
nTime = 193 #total number of movie TRs
boundaries <- read.csv('../../behavioral-data/CamCan_movie_event_onsets.csv')  # from Ben-Yakov & Henson (2018)
boundaries$TR <- boundaries$Boundary/my.TR

# I am now classifying time points as in (1) or outside (0) of an event transition window,
# with a 5-TR window centered on each boundary
boundary.windows <- array(0,c(nTime))
boundaries$adjustedTR <- round(boundaries$TR + 2)  #accounting for HRF lag (~5s) before rounding to nearest TR
for (b in 1:nrow(boundaries)){
  min.point <- boundaries$adjustedTR[b] - 2
  max.point <- boundaries$adjustedTR[b] + 2
  boundary.windows[min.point:max.point] <- 1
}


### option to re-define nTime if we only want to analyze up to a certain point. 
### added to check the influence of the "gun shot" event toward the end of the movie (no change to results)  
# nTime = 175
boundary.windows <- boundary.windows[1:nTime]

cat('Analyzing PMN connectivity over',nTime,'TRs (',round(nTime*my.TR,2),'seconds ) of movie watching\n')
## Analyzing PMN connectivity over 193 TRs ( 476.71 seconds ) of movie watching

Data Inspection

ROI time series

Mean (across voxels) time series from unsmoothed data

timeseries.file <- paste('./',this.group,'/PM_node_timeseries_',this.group,'.csv',sep="")
cat('Loading PMN time series from:',timeseries.file,'\n')
## Loading PMN time series from: ./Discovery/PM_node_timeseries_Discovery.csv
allData = read.csv(timeseries.file, header = TRUE)
# for testing -- analyze subset of time points
allData <- subset(allData, Time <=nTime)


allData$Subject=as.factor(allData$Subject)
allData$Node=as.factor(allData$Node)

subjects <- levels(allData$Subject)
NSubs    <- length(subjects)
cat('Total Number of Subjects =',NSubs,'\n')
## Total Number of Subjects = 68
rois = levels(allData$Node)
cat('ROIs =',rois,'\n')
## ROIs = aAG Hipp MPFC pAG PCC PHC Prec RSC
# number of unique connections for 8x8 rois
nConnections <- ((length(rois)*length(rois)) - length(rois))/2

To facilitate data inspection, here is an interactive plot showing the denoised, mean time series of each ROI per subject.
Each time series has been standardized within subject and ROI.

scaled_data <- allData %>% group_by(Subject, Node) %>%
                   mutate(Z = scale(Value))

lim <- ceil(max(abs(scaled_data$Z)))

plot_ly(scaled_data, x = ~Time, y = ~Z, frame = ~Subject, color = ~Node,
        colors = mycolors, type = "scatter", mode = "lines",
        width = 900, height = 400) %>%
  layout(title = "Z scored PM time series",
         xaxis = list(title = "TR"),
         yaxis = list(range = c(-lim,lim)))  %>%
  animation_opts(frame = 500, transition = 100, redraw = FALSE)

Mean standardized time series across subjects, by ROI:

summary_timeseries <- scaled_data %>% group_by(Time, Node) %>%
                         summarise(MeanZ = mean(Z),
                                   SEZ = se(Z))

ggplot(summary_timeseries) +
  scale_fill_manual(values = mycolors) +
  scale_color_manual(values = mycolors) +
  geom_hline(yintercept=0, size=1) +
  geom_ribbon(alpha=0.3, aes(x=Time*my.TR, fill=Node,
                             ymin=MeanZ-SEZ, ymax=MeanZ+SEZ), color=NA) +
  geom_line(aes(x=Time*my.TR, y=MeanZ, color=Node), size=1.5) +
  scale_y_continuous(expand=c(0,0.01,0,0), limits=c(-2,2.5)) +
  scale_x_continuous(expand=c(0,0,0,0)) +
  labs(x="Time (seconds)", y='Mean Z') +
  ggtitle('Group-averged ROI time-series') +
  theme(plot.title = element_text(size = 30, face='plain', hjust=0.5),
        axis.line.x = element_line(color = "black", size = 1),
        axis.line.y = element_line(color = "black", size = 1),
        axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 18),
        axis.title.y = element_text(size = 22), axis.title.x = element_text(size = 22),
        panel.background = element_blank(),
        legend.title = element_blank(),
        legend.key.size = unit(2,"line"),
        legend.margin = margin(0,0,0,0, "cm"),
        plot.margin = margin(0.25, 0.25, 0.25, 0.25, "cm"))

Framewise displacement

Even though FD, and other confounds, have been removed from the above ROI time series, here’s the TR-by-TR FD for reference:

fdData = read.csv('../../dataset-info/FD_byTR.csv', header = TRUE)

# filter by subjects in this sample:
fdData <- fdData[fdData$SubID %in% subjects,] %>% 
               droplevels()
fdData$SubID=as.factor(fdData$SubID)
colnames(fdData)[2] <- 'Time'

fdData <- subset(fdData, Time <= nTime)


summary_FD <- fdData %>% group_by(Time) %>%
                         summarise(Mean = mean(FD),
                                   SE = se(FD))

ggplot(summary_FD) +
  geom_ribbon(alpha=0.3, aes(x=Time*my.TR, ymin=Mean-SE, ymax=Mean+SE), 
              color=NA, fill="gray40") +
  geom_line(aes(x=Time*my.TR, y=Mean), size=1.5, color="gray20") +
  scale_y_continuous(expand=c(0,0,0,0), limits=c(0,.24)) +
  scale_x_continuous(expand=c(0,0,0,0)) +
  labs(x="Time (seconds)", y='Mean FD') +
  ggtitle('Group-averged FD') +
  theme(plot.title = element_text(size = 30, face='plain', hjust=0.5),
        axis.line.x = element_line(color = "black", size = 1),
        axis.line.y = element_line(color = "black", size = 1),
        axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 18),
        axis.title.y = element_text(size = 22), axis.title.x = element_text(size = 22),
        panel.background = element_blank(),
        legend.title = element_blank(),
        legend.key.size = unit(2,"line"),
        legend.margin = margin(0,0,0,0, "cm"),
        plot.margin = margin(0.25, 0.25, 0.25, 0.25, "cm"))

Seed-to-voxel connectivity

Whole-brain data was smoothed (6mm FWHM) and masked with the MNI brainmask from normalization. Calculates the top 20% (range 30-10%) of voxel connections for each ROI.
Note. Whole-brain connectivity was calculated in MATLAB, and is loaded in here as a csv file and pre-generated BrainNet Viewer images.

th_WB = 0.80 #top 20 % of connections -- for visualization
# also iterate over other roi-specific network densities to check the stability of communities
wb_iter = seq(from=.7, to=.9, by=.05) # 30% to 10% network densitiies


# read in csv with group-averaged seed-to-voxel connectivity patterns (r values)
seedtovoxel.file <- paste('./',this.group,'/wholebrain/group-average/group_PM_wholebrain_connectivity.csv',sep="")
cat('Loading seed-to-voxel connectivity values from:',seedtovoxel.file,'\n')
## Loading seed-to-voxel connectivity values from: ./Discovery/wholebrain/group-average/group_PM_wholebrain_connectivity.csv
roi_wholebrain   <- read.csv(seedtovoxel.file)


# order columns by my roi order (currently alphabetical):
ord = match(rois,colnames(roi_wholebrain))
roi_wholebrain <- roi_wholebrain[,ord]

# store binarized whole brain connections (from group-level averaged connectivity) per seed, per density threshold
wholebrain_th <- array(NA,c(nrow(roi_wholebrain),ncol(roi_wholebrain),length(wb_iter)))
colnames(wholebrain_th) <- rois
for (d in 1:length(wb_iter)) {
  for (r in 1:length(rois)) {
    th <- quantile(roi_wholebrain[,r],wb_iter[d])
    wholebrain_th[roi_wholebrain[,r] >= th,r,d] <- 1
    wholebrain_th[roi_wholebrain[,r] < th,r,d]  <- 0
  }
}

# plot pre-generated images of each roi's top 20% of connections:  
plots.images = list()
for (r in 1:length(rois)) {
  plots.images[[r]] <- ggdraw() +
                       draw_image(paste('./',this.group,'/wholebrain/group-average/task-movie_',rois[r],
                                        '_R_seedtovoxel_binarized_quantile_',th_WB,'.png',sep="")) + 
                       ggtitle(rois[r]) +
                       theme(plot.title = element_text(hjust = 0.5, size=20, margin=margin(0,0,5,0)),
                             plot.margin = margin(0, 0, 0, 0, "cm"))
}

Now, I am using the Louvain method to detect ROI clusters based on the similarity of whole-brain connectivity patterns.

  • Default gamma = 1 – using a range of values and then taking the % of iterations that each pair of ROIs is placed into the same module.
  • Group with hierarchical clustering (hclust).
gammas = seq(0.75,1.25,0.01) #default Louvain gamma = 1

group_modules <- data.frame(array(0,c(length(gammas)*length(wb_iter),length(rois))))
colnames(group_modules) <- rois

n = 0
# Calculate Louvain modules over densities...
for (d in 1:length(wb_iter)) {
  # correlate wholebrain seed-to-voxel connectivity patterns between rois
  net <- cor(wholebrain_th[,,d])

  # ... and gammas
  for (g in gammas) {
    n = n+1
    modules <- louvain(net,g)$community
    group_modules[n,] <- modules
  }
} # end of loop over densities


# calculate % of the time each ROI is grouped with every other (across all density + gamma levels):
roi_modules <- array(0,c(length(rois),length(rois)))
colnames(roi_modules) <- rois
rownames(roi_modules) <- rois
for (x in 1:length(rois)) {
  x_vec <- group_modules[,colnames(group_modules) == rois[x]]
  for (y in 1:length(rois)) {
    y_vec <- group_modules[,colnames(group_modules) == rois[y]]
    xy <- x_vec == y_vec
    roi_modules[x,y] <- (sum(xy)/length(x_vec))*100
  }
}



### plot ------
my.modules <- as.vector(group_modules[30,]) # manually selecting the best parcellation after visual inspection of clustering
n.mod <- max(unique(my.modules))



########## difference between discovery and replication markdowns #########
if (this.group == "Discovery") {
  # sort by communities
  ord = c()
  for (n in 1:n.mod) {
    idx <- which(my.modules == n)
    
    #hclust within module in case any more fine-grained module groupings
    h <- hclust(as.dist(100-roi_modules[idx,idx]))
    ord <- c(ord, idx[h$order])
  }
  
  # save this roi order for use with replication ample
  save(ord, file='roi_order.RData')
  
} else if (this.group == "Replication") {
  # load in ROI order from discovery sample to keep visualization consistent
  load('roi_order.RData')  
}
############################################################################



# ***** re-order by community ******* #
roi_modules  <- roi_modules[ord,ord]
my.modules   <- my.modules[ord]
allData$Node <- factor(allData$Node,levels(allData$Node)[ord])
mycolors     <- mycolors[ord]
rois         <- rois[ord]
# *********************************** #


# plot module matrix
data <- melt(roi_modules)
# remove diagonal
data <- data[data$Var1 != data$Var2,]

p2 <- ggplot(data, aes(Var1, Var2, fill=value)) +
  geom_tile(color="black", size=.5) +
  geom_text(aes(label = round(value)), size = 4, color="white") +
  scale_fill_gradientn(limits=c(40,100), colors = gray.colors(12), na.value = "white",
                       guide = guide_colorbar(frame.colour = "black")) +
  labs(fill="% Same\nModule", tag="c") +
  ggtitle('Community probability')+
  theme(plot.title = element_text(hjust = 0.5, size=24, margin=margin(0,0,20,0), face="plain"),  
        axis.text.x = element_text(size = 16, color = mycolors, hjust=1, vjust=1, angle=45), 
        axis.text.y = element_text(size = 16, color = mycolors),
        axis.ticks.x = element_blank(),axis.ticks.y = element_blank(),
        axis.line.x = element_blank(),axis.line.y = element_blank(),
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        legend.title = element_text(size=14), legend.text = element_text(size=12),
        legend.key.width = unit(0.4,"cm"), legend.key.height = unit(0.8,"cm"),
        plot.margin = margin(0, 0, 0, 0, "cm"),
        plot.tag = element_text(size = 34, face = "bold"))

# add rectangles to highlight communities
min.x <- .5
for (n in 1:n.mod) {
  count.n <- sum(my.modules == n)
  max.x <- min.x + count.n
  p2 <- p2 +
    geom_rect(xmin=min.x, xmax=max.x, ymin=min.x, ymax=max.x,
              fill=NA, color=module.colors[n],
              size=2)
  min.x <- max.x
}



# add images showing overlap of connections by submodule
p3 <- ggdraw() +
  draw_image(paste('./',this.group,'/wholebrain/group-average/task-movie_quantile_',th_WB,'_sum_Core.png',sep="")) + 
  ggtitle('"Core": 20% density') + labs(tag="d") +
  geom_label(aes(x=0.5, y=0.065, label="Connected Regions"), 
             size=4.5, color="black", fill="white", label.size=0) +
  geom_label(aes(x=0.5, y=0.115, label="0                        5"), 
             size=4, color="black", fill=NA, label.size=0) +
  theme(plot.title = element_text(hjust = 0.5, size=24, margin=margin(0,0,10,0)),
        plot.tag = element_text(size = 34, face="bold"),
        plot.margin = margin(0, 0, 0, 0, "cm"))

p4 <- ggdraw() +
  draw_image(paste('./',this.group,'/wholebrain/group-average/task-movie_quantile_',th_WB,'_sum_Ventral.png',sep="")) + 
  ggtitle('"Ventral": 20% density') + labs(tag="e") +
  geom_label(aes(x=0.5, y=0.065, label="Connected Regions"), 
             size=4.5, color="black", fill="white", label.size=0) +
  geom_label(aes(x=0.5, y=0.115, label="0                        3"), 
             size=4, color="black", fill=NA, label.size=0) +
  theme(plot.title = element_text(hjust = 0.5, size=24, margin=margin(0,0,10,0)),
        plot.tag = element_text(size = 34, face="bold"),
        plot.margin = margin(0, 0, 0, 0, "cm"))
# add in roi figure and analysis schematic ***** only for Discovery sample ******
analysis.plot <- ggdraw() +
  draw_image("./methods-figs/method_1_fig.png") + 
  labs(tag="a") +
  theme(plot.tag = element_text(size = 34, face="bold"),
        plot.margin = margin(0, 0, 0, 0, "cm"))


# re-order panel A by module assignments and group:
p1 <- ggarrange(plotlist=plots.images[ord], ncol = 4, nrow = 2) +
          labs(tag="b") + theme(plot.tag = element_text(size = 34, face="bold"),
                                plot.margin = margin(0, 0, 0, 0, "cm"))

top.plot <- ggarrange(analysis.plot,p1, ncol = 2, nrow = 1, widths=c(1,1))


###### combine whole brain plots
p <- ggarrange(p2,p3,p4, ncol = 3, nrow = 1, widths=c(1.05,1,1))
myPlot <- ggarrange(top.plot,p, ncol = 1, nrow = 2, heights=c(1.25,1))
plot(myPlot)

## save
ggsave(paste('./',this.group,'/figures/seed-to-voxel_',this.group,'.pdf',sep=""), plot = myPlot, device = "pdf", width=15, height=10.5, unit="in")

From now on, ROIs are sorted into the above modules – ‘Ventral’ and ‘Core’.

Intranetwork connectivity

Bivariate network = Pearson’s correlation between the time-series of PMM regions, per subject.
Partial network = correlation between ROIi and ROIj and controlling for all 6 other nodes.

## create full pearson's correlation connectivity matrix --> pearson's r
connMatrix = array(data = 0, c(length(rois),length(rois),NSubs))
colnames(connMatrix) <- rois
rownames(connMatrix) <- rois
connMatrix.partial   <- connMatrix

for (s in 1:NSubs) {
  
  subData   <- subset(allData, Subject == subjects[s])
  # from long to wide format
  subMatrix <- spread(subData, Node, Value)
  connMatrix[,,s]         <- cor(subMatrix[,match(rois,colnames(subMatrix))])
  connMatrix.partial[,,s] <- cor2pcor(connMatrix[,,s]) #convert bivariate to partial network
  diag(connMatrix[,,s]) <- 0
  diag(connMatrix.partial[,,s]) <- 0

} #end of loop through subjects


# store
node.bivariate <- connMatrix
node.partial   <- connMatrix.partial


### calculate mean for each connection across subjects, applying fisher z transformation before averaging, 
### with group-values converted back to r
node.bivariate.mean  <- fisherz2r(apply(fisherz(node.bivariate), c(1,2), mean))
node.partial.mean    <- fisherz2r(apply(fisherz(node.partial), c(1,2), mean))

Average connectivity within and between PMN subsystems:

  • Bivariate:
# labels connections by module:
data.mod <- melt(node.bivariate)  #Var3 = Subject index
data.mod <- data.mod[data.mod$value != 0,] # remove diagonal
data.mod$Module = ''
data.mod$Module[data.mod$Var1 %in% rois[my.modules==1] & data.mod$Var2 %in% rois[my.modules==1]] <- 'Core'
data.mod$Module[data.mod$Var1 %in% rois[my.modules==2] & data.mod$Var2 %in% rois[my.modules==2]] <- 'Ventral'
data.mod$Module[(data.mod$Var1 %in% rois[my.modules==1] & data.mod$Var2 %in% rois[my.modules==2]) |
                (data.mod$Var1 %in% rois[my.modules==2] & data.mod$Var2 %in% rois[my.modules==1])] <- 'Ventral-Core'


# within-subject mean module connectivity (fisher z transformed)
data.summary.subject <- data.frame(
                          data.mod %>% 
                          group_by(Var3, Module) %>%
                          summarise(Strength = mean(fisherz(value)))
                          )


# plot subject-level averages
p1 <- ggplot(data.summary.subject, aes(x=Module, y=fisherz2r(Strength), color=Module)) +
  scale_fill_manual(values = module.colors) +
  scale_color_manual(values = module.colors) +
  geom_point(size = 2, position=position_jitter(width = .15), alpha=.6, color="gray60") +
  geom_boxplot(size=1, width=.6, fill=NA, outlier.color=NA) +
  geom_hline(yintercept = 0, size = 1) +
  scale_y_continuous(expand = c(0,0.01,0,0.01), limits=c(-.25,.65)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 7)) +
  labs(y='R', tag="c") +
  theme(axis.line.x = element_line(color = "black", size = 1), axis.line.y = element_line(color = "black", size = 1),
        axis.text.x = element_text(size = 20, angle=30, hjust=1, vjust=1),
        axis.text.y = element_text(size = 18),
        axis.title.y = element_text(size = 22), axis.title.x = element_blank(),
        panel.background = element_blank(),
        legend.position="none",
        plot.tag = element_text(size = 34, face="bold"),
        plot.margin = margin(0.1, 0.1, 0, 0.1, "cm"))


#### means and CIs per module
data.summary.group <- data.summary.subject %>% group_by(Module) %>%
                           summarise(Mean = mean(Strength),
                                     CI = ci(Strength))
kable(data.summary.group, format = "pandoc",
      caption = "Mean bivariate connectivity within and between PM modules")
Mean bivariate connectivity within and between PM modules
Module Mean CI
Core 0.3453571 0.0240332
Ventral 0.3982606 0.0340449
Ventral-Core 0.2889817 0.0219783 
#### t-tests to validate whole-brain defined modules in terms of bivariate connectivity
# runs on each pair of modules
ts <- data.summary.subject %>%
            pairwise_t_test(
                Strength ~ Module, paired = TRUE, 
                p.adjust.method = "none",
                detailed=T
            )
kable(ts, format = "pandoc",
      caption = "Bivariate connectivity differences within and between PM modules")
Bivariate connectivity differences within and between PM modules
estimate .y. group1 group2 n1 n2 statistic p df conf.low conf.high method alternative p.adj p.adj.signif
-0.0529035 Strength Core Ventral 68 68 -2.853613 0.0060000 67 -0.0899078 -0.0158992 T-test two.sided 0.0060000 **
0.0563754 Strength Core Ventral-Core 68 68 4.792147 0.0000095 67 0.0328941 0.0798567 T-test two.sided 0.0000095 ****
0.1092789 Strength Ventral Ventral-Core 68 68 6.630783 0.0000000 67 0.0763836 0.1421742 T-test two.sided 0.0000000 **** 
### save out time-averaged module connectivity for episodic memory individual differences analyses
data.summary.subject$Var3 <- subjects[data.summary.subject$Var3] # add in actual subject IDs from indexes
timeaveraged.file <- paste('./',this.group,'/module-connectivity_',this.group,'.RData',sep="")
cat('Saving time-averaged module connectivity to:',timeaveraged.file,'\n')
## Saving time-averaged module connectivity to: ./Discovery/module-connectivity_Discovery.RData
save(data.summary.subject, file=timeaveraged.file)
  • Partial:
# labels connections by module:
data.mod <- melt(node.partial)  #Var3 = Subject index
data.mod <- data.mod[data.mod$value != 0,] # remove diagonal
data.mod$Module = ''
data.mod$Module[data.mod$Var1 %in% rois[my.modules==1] & data.mod$Var2 %in% rois[my.modules==1]] <- 'Core'
data.mod$Module[data.mod$Var1 %in% rois[my.modules==2] & data.mod$Var2 %in% rois[my.modules==2]] <- 'Ventral'
data.mod$Module[(data.mod$Var1 %in% rois[my.modules==1] & data.mod$Var2 %in% rois[my.modules==2]) |
                (data.mod$Var1 %in% rois[my.modules==2] & data.mod$Var2 %in% rois[my.modules==1])] <- 'Ventral-Core'


# within-subject mean module connectivity (fisher z transformed)
data.summary.subject <- data.frame(
                          data.mod %>% 
                          group_by(Var3, Module) %>%
                          summarise(Strength = mean(fisherz(value)))
                          )


# plot subject-level averages
p2 <- ggplot(data.summary.subject, aes(x=Module, y=fisherz2r(Strength), color=Module)) +
  scale_fill_manual(values = module.colors) +
  scale_color_manual(values = module.colors) +
  geom_point(size = 2, position=position_jitter(width = .15), alpha=.6, color="gray60") +
  geom_boxplot(size=1, width=.6, fill=NA, outlier.color=NA) +
  geom_hline(yintercept = 0, size = 1) +
  scale_y_continuous(expand = c(0,0.01,0,0.01), limits=c(-.25,.65)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 7)) +
  labs(y='R', tag="d") +
  theme(axis.line.x = element_line(color = "black", size = 1), axis.line.y = element_line(color = "black", size = 1),
        axis.text.x = element_text(size = 20, angle=30, hjust=1, vjust=1),
        axis.text.y = element_text(size = 18),
        axis.title.y = element_text(size = 22), axis.title.x = element_blank(),
        panel.background = element_blank(),
        legend.position="none",
        plot.tag = element_text(size = 34, face="bold"),
        plot.margin = margin(0.1, 0.1, 0, 0.1, "cm"))


#### means and CIs per module
data.summary.group <- data.summary.subject %>% group_by(Module) %>%
                           summarise(Mean = mean(Strength),
                                     CI = ci(Strength))
kable(data.summary.group, format = "pandoc",
      caption = "Mean partial connectivity within and between PM modules")
Mean partial connectivity within and between PM modules
Module Mean CI
Core 0.1308432 0.0076399
Ventral 0.2153361 0.0232042
Ventral-Core 0.0721908 0.0080519 
#### t-tests to validate whole-brain defined modules in terms of bivariate connectivity
# runs on each pair of modules
ts <- data.summary.subject %>%
            pairwise_t_test(
                Strength ~ Module, paired = TRUE, 
                p.adjust.method = "none",
                detailed=T
            )
kable(ts, format = "pandoc",
      caption = "Partial connectivity differences within and between PM modules")
Partial connectivity differences within and between PM modules
estimate .y. group1 group2 n1 n2 statistic p df conf.low conf.high method alternative p.adj p.adj.signif
-0.0844929 Strength Core Ventral 68 68 -7.590603 0 67 -0.1067109 -0.0622748 T-test two.sided 0 ****
0.0586524 Strength Core Ventral-Core 68 68 8.024649 0 67 0.0440635 0.0732412 T-test two.sided 0 ****
0.1431452 Strength Ventral Ventral-Core 68 68 9.530374 0 67 0.1131654 0.1731251 T-test two.sided 0 ****

Mean bivariate and partial networks as graphs:

scale.factor <- 8  #just to scale edge weights for visualization, no impact on stats
c.th = .2  # threshold for visualization of 'strong' connections - based on our time points, it reflects critical r at p < .005
p.th = 0   # testing for significance against 0
###########


# bivariate --------------------------------------------------------------------
#plot those that are p-corrected greater than 0
ps <- apply(fisherz(node.bivariate), c(1,2),
            function(x) t.test(x, mu = p.th, alternative="greater")$p.value)
ps <- ps * nConnections # bonferroni

net <- node.bivariate.mean
net[ps >= .05] = 0  #remove non-significant
net[net < 0]   = 0  #remove negative (should be redundant with above)


#plot network, highlighting strong connections
net <- net*scale.factor
net <- network(net, directed = FALSE, ignore.eval = FALSE, names.eval = "weights")
network::set.edge.attribute(net, "color",
                            ifelse(net %e% "weights" > c.th*scale.factor, "gray50", "gray80"))
network.vertex.names(net) = rois
p3 <- suppressMessages(ggnet2(net,
                 node.size = 11, node.color = mycolors,
                 label=TRUE, label.size=5,
                 edge.size = "weights", edge.color = "color",
                 layout.exp = 0.05) +
  scale_y_continuous(expand=c(.05,.05,.05,.05)) +
  scale_x_continuous(expand=c(.05,.05,.05,.05)) +
  theme(plot.margin = margin(0.1, 0.1, 0.1, 1.1, "cm"),
        axis.line.x = element_blank(), axis.line.y = element_blank(),
        axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        axis.text.x = element_blank(), axis.ticks.x = element_blank()))

p3 <- ggarrange(p3) + 
      labs(tag="a") + 
      theme(plot.tag = element_text(size = 34, face="bold"),
            plot.margin = margin(0.1, 0, 0, 0.1, "cm"))



#partial ----------------------------------------------------------------------
#plot those that are p-corrected greater than 0
ps <- apply(fisherz(node.partial), c(1,2),
            function(x) t.test(x, mu = p.th, alternative="greater")$p.value)
ps <- ps * nConnections # bonferroni

net <- node.partial.mean
net[ps >= .05] = 0  #remove non-significant
net[net < 0]   = 0  #remove negative (should be redundant with above)


#plot network, highlighting strong connections
net <- net*scale.factor
net <- network(net, directed = FALSE, ignore.eval = FALSE, names.eval = "weights")
network::set.edge.attribute(net, "color",
                            ifelse(net %e% "weights" > c.th*scale.factor, "gray50", "gray80"))
network.vertex.names(net) = rois
p4 <- suppressMessages(ggnet2(net,
                 node.size = 11, node.color = mycolors,
                 label=TRUE, label.size=5,
                 edge.size = "weights", edge.color = "color",
                 layout.exp = 0.05) +
  scale_y_continuous(expand=c(.05,.05,.05,.05)) +
  scale_x_continuous(expand=c(.05,.05,.05,.05)) +
  theme(plot.margin = margin(0.1, 0.1, 0.1, 1.1, "cm"),
        axis.line.x = element_blank(), axis.line.y = element_blank(),
        axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        axis.text.x = element_blank(), axis.ticks.x = element_blank()))

p4 <- ggarrange(p4) + 
      labs(tag="b") + 
      theme(plot.tag = element_text(size = 34, face="bold"),
            plot.margin = margin(0.1, 0, 0, 0.1, "cm"))
myPlot <- ggarrange(p3,p4,p1,p2, ncol = 2, nrow=2, heights=c(1,1.1))
plot(myPlot)

## save as Rdata for merging with other group
save(myPlot, file=paste('./',this.group,'/figures/intranetwork_',this.group,'.RData',sep=""))

Virtual lesion

Here, I’m testing the mediating influence of each node on other connections within the network – if we partial out (remove) each node in turn rather than all at once (partial network, above), how do the other connections in the network change?

For each node, I’m generating a circular graph without that node and it’s connections. The width of each edge represents the strength of the original connection (c in a mediation model). The color of the edge (from gray to blue) represents how much that connection has decreased after ‘lesioning’ a node (Pm = 1-[c’/c], where c’ is the partial correlation in a mediation model). [Pm = proportion mediated]

### group-level bivariate mask -- 
### only plot connections originally > threshold for visualization
### (to make sure there is something meaningful to mediate)
b.mask <- node.bivariate.mean
b.mask[b.mask < c.th] <- 0
b.mask[b.mask > 0] <- 1
diag(b.mask) <- 0

# color gradient of edge colors to show % connection strength mediated
edge.bins    <- seq(.1,1,by=.1) # increments of Pm
edge.colors  <- blue.colors(length(edge.bins))
# ---------------------------------------------------------------------# 


# to store summary network Pm stats:
roi_effect = data.frame(array(0,c(length(rois)*NSubs,3)))
colnames(roi_effect) = c('SubID','Node','Influence')
row = 0

# to store partial networks after removing RSC and PCC, to illustrate strongest effects
lesion.examples <- array(data = NA, c(length(rois),length(rois),2))
ex = 0


plots.graphs = list()
for (p in 1:length(rois)) {
  #for full bivariate (ignoring this node's connections)
  connMatrix.c = array(data = NA, c(length(rois),length(rois),NSubs))
  colnames(connMatrix.c) <- rois
  rownames(connMatrix.c) <- rois
  #for partial (removing mediating influence of this roi)
  connMatrix.p <- connMatrix.c  

  for (s in 1:length(subjects)) {
    subData <- subset(allData, Subject == subjects[s])
    #time-series for this node (mediator)
    pData   <- subData$Value[subData$Node == rois[p]] 
    
    for (y in 1:length(rois)) {
      #time-series for ROIy
      yData = subData$Value[subData$Node == rois[y]]
      for (x in 1:length(rois)) {
        # time-series for ROIx
        xData = subData$Value[subData$Node == rois[x]]
        if ((y != p) & (x != p) & (x != y)) {
          # partial correlation between x and y when controlling for p
          connMatrix.p[y,x,s] <- pcor(cbind(xData,yData,pData))$estimate["xData","yData"]
          # original, bivariate correlation between x and y
          connMatrix.c[y,x,s] <- cor(xData,yData)
        }
      }
    }
    
    
    # per subject, what is the % of variance in all other connections accounted for? 
    # (% change from mean original PMN correlation to mean partial)
    row = row + 1
    roi_effect$SubID[row] <- s  #subject index
    roi_effect$Node[row]  <- rois[p]
    p.conn <- fisherz2r(mean(fisherz(connMatrix.p[,,s]), na.rm = TRUE)) # average partial correlation (c')
    c.conn <- fisherz2r(mean(fisherz(connMatrix.c[,,s]), na.rm = TRUE)) # average original correlation (c)
    Pm <- 1 - (p.conn/c.conn) # proportion mediated
    if (Pm < 0) { #mask any small "supressing" effects where average p.conn is actually larger than average c.conn
      Pm <- 0
    } else if (Pm > 1) { #if the average p.conn is now negative
      Pm <- 1
    }
    roi_effect$Influence[row] <- Pm

  } #end of loop through subjects


  ### calculate group-level Pm for each PMN connection ----------------------- #
  connmean.p <- fisherz2r(apply(fisherz(connMatrix.p), c(1,2), mean))
  ## store partial network for example plot if PCC or RSC
  if (rois[p] == "PCC" || rois[p] == "RSC") {
    ex = ex + 1
    lesion.examples[,,ex] <- connmean.p
  }
  connmean.p[b.mask == 0] <- NA #ignore connections that weren't present above our threshold (needs to be something to reduce)
  connmean.c <- fisherz2r(apply(fisherz(connMatrix.c), c(1,2), mean))
  connmean.c[b.mask == 0] <- NA
  
  # pm for each connection
  connmean.r <- 1 - (connmean.p/connmean.c)
  connmean.r[connmean.r < 0] <- 0.001 #mask for edges where there is a small increase for partial cor ("supressing" effects)
  # ^^ setting non-zero just so it gets recognized as an edge below
  connmean.r[connmean.r > 1] <- 1     #where p.conn is now negative, all variance in c.conn is explained. 

  

  ##### PLOT ---------------------------------------------------------------- #

  #specify initial network of variance explained (Pm) to get edge colors
  net <- connmean.r
  net[is.na(net)] <- 0
  net <- network(net, directed = FALSE, ignore.eval = FALSE, names.eval = "weights")

  # set gradient color for % reduction per edge
  my.edges     <- network::get.edge.value(net, "weights")
  edge.colors.sort <- array('',c(length(my.edges)))
  for (e in 1:length(my.edges)) {
    diff <- edge.bins - my.edges[e]
    diff[diff < 0] = 2  #set above max to get next lowest bin
    edge.colors.sort[e] <- edge.colors[diff == min(diff)]
  }
  # if it's the first mediator ROI, get a color scale to later add to the plot:
  if (p == 1) {
    net.df <- ggnetwork(net, layout="circle")
    edge.p <- ggplot(net.df, aes(x = x, y = y, xend = xend, yend = yend,
                            color=weights)) +
      geom_edges() +
      scale_color_gradientn(limits = c(min(edge.bins),max(edge.bins)),
                            colors  = edge.colors, na.value = "white",
                            guide = guide_colorbar(frame.colour = "black")) +
      labs(color="Proportion\nmediated") +
      theme(legend.text = element_text(size=16), 
            legend.title = element_text(size=18, margin = margin(0,0,5,0)),
            legend.key.width = unit(0.4,"cm"), legend.key.height = unit(1,"cm"),
            plot.margin = margin(0, 0.25, 0, 0, "cm"))
    edge.legend <- as_ggplot(get_legend(edge.p))
  }
  
  
  #now replace network with the original bivariate connections (for edge width)
  net <- connmean.c*6  #scale just for visualization purposes
  net[is.na(net)] <- 0
  net <- network(net, directed = FALSE, ignore.eval = FALSE, names.eval = "weights")
  
  # draw graph
  # make current node white so not visible (removed statistically)
  lesion.colors <- mycolors
  lesion.colors[p] <- "white" #remove current mediator node
  network.vertex.names(net) = rois
  plots.graphs[[p]] <- suppressMessages(ggnet2(net, mode = "circle",
                              node.size = 10, node.color = lesion.colors,
                              label=FALSE,
                              edge.size = "weights", edge.color = edge.colors.sort,
                              layout.exp = 0.05) +
    scale_y_continuous(expand=c(.05,.05,.05,.05)) +
    scale_x_continuous(expand=c(.05,.05,.05,.05)) +
    theme(plot.margin = margin(0.25, 0.5, 1, 0.5, "cm"),
          plot.title = element_text(size = 28, margin=margin(5,0,10,0), face="plain", hjust=0),
          axis.line.x = element_blank(), axis.line.y = element_blank(),
          axis.text.y = element_blank(), axis.ticks.y = element_blank(),
          axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
    ggtitle(rois[p]) + border(size = 1, linetype = 1))


}  # end of loop through mediator rois ---------------------------------------------- #

Proportion of network connectivity mediated, by ROI:

## Add summary node influence bar graph ---------------------------------------
roi_effect$Node <- as.factor(roi_effect$Node)
roi_effect$Node <- factor(roi_effect$Node,levels(roi_effect$Node)[ord])
roi_effect.summary <- roi_effect %>%
                       group_by(Node) %>%
                       summarise(Strength = mean(Influence), CI = ci(Influence))

# now reorder rois (and their colors) by mean mediating effect for these plots
med.ord    <- order(roi_effect.summary$Strength)
med.rois   <- rois[med.ord]
med.colors <- mycolors[med.ord]
roi_effect$Node <- factor(roi_effect$Node,levels(roi_effect$Node)[med.ord])


### plot average influence per roi:
p1 <- ggplot(roi_effect, aes(x=Node, y=Influence, color=Node, fill=Node)) +
     scale_fill_manual(values = med.colors) +
     scale_color_manual(values = med.colors) +
     geom_point(size = 2, position=position_jitter(width = .1), alpha=.6, color="gray60") +
     geom_boxplot(size=1, width=.6, fill=NA, outlier.color=NA) +
     scale_y_continuous(expand = c(0,0,0,0), limits=c(-.01,1.01)) +
     labs(y='Network Pm', tag="c") +
     theme(axis.line.x = element_line(color = "black", size = 1),
          axis.line.y = element_line(color = "black", size = 1),
          axis.text.x = element_text(size = 22, hjust=1, vjust=1, angle=45),
          axis.text.y = element_text(size = 20),
          axis.title.y = element_text(size = 22), axis.title.x = element_blank(),
          panel.background = element_blank(),
          legend.position="none",
          plot.tag = element_text(size = 40, face="bold"),
          plot.margin = margin(0.25, 0.5, 0.25, 0, "cm"))

Illustrate the resulting network (thresholded at r > .2) if you remove RSC or PCC (largest Pm):

lesion.examples[is.na(lesion.examples)] <- 0
lesion.examples[lesion.examples < c.th] <- 0  #remove low correlations


### RSC
net <- lesion.examples[rois!="RSC",rois!="RSC",2]
net <- network(net, directed = FALSE, ignore.eval = FALSE, names.eval = "weights")

# draw graph
network.vertex.names(net) = rois[rois!="RSC"]
p.lesionA <- suppressMessages(ggnet2(net,
                    node.size = 8, node.color = mycolors[rois!="RSC"],
                    label=FALSE,
                    edge.size = 2, edge.color = "gray70",
                    layout.exp = 0.1) +
  scale_y_continuous(expand=c(.05,.05,.05,.05)) +
  theme(plot.margin = margin(0, 0.5, 0, 3, "cm"),
        axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  border(size = 2, linetype = 1, color=mycolors[rois=="RSC"]))


### PCC
net <- lesion.examples[rois!="PCC",rois!="PCC",1]
net <- network(net, directed = FALSE, ignore.eval = FALSE, names.eval = "weights")

# draw graph
network.vertex.names(net) = rois[rois!="PCC"]
p.lesionB <- suppressMessages(ggnet2(net,
                    node.size = 8, node.color = mycolors[rois!="PCC"],
                    label=FALSE,
                    edge.size = 2, edge.color = "gray70",
                    layout.exp = 0.1) +
  scale_y_continuous(expand=c(.05,.05,.05,.05)) +
  theme(plot.margin = margin(0, 1, 0, 0.5, "cm"),
        axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  border(size = 2, linetype = 1, color=mycolors[rois=="PCC"]))


# join
p.lesion <- ggarrange(p.lesionA,p.lesionB, ncol=2, widths=c(1.32,1)) +
            labs(tag="b") + theme(plot.tag = element_text(size = 40, face="bold"))
# order and plot graph panels:
myPlot <- ggarrange(plotlist=plots.graphs[med.ord], ncol = 4, nrow = 2) +
          labs(tag="a") + theme(plot.tag = element_text(size = 40, face="bold"))

# add lesion examples over summary plot
myPlotB <- ggarrange(p.lesion, p1, nrow=2, heights=c(2,5))

# join! (with edge color legend)
joined.Plot <- ggarrange(myPlot, edge.legend, myPlotB, ncol=3, nrow=1, widths=c(2,0.25,1))
plot(joined.Plot)

## save as Rdata for merging with other group
save(joined.Plot, file=paste('./',this.group,'/figures/virtual-lesion_',this.group,'.RData',sep=""))